# !/usr/bin/env python
# -*- coding: utf-8 -*-
"""
@Time        : 2021/1/3 11:23
@Author      : Albert Darren
@Contact     : 2563491540@qq.com
@File        : numerical_integral.py
@Version     : Version 1.0.0
@Description : TODO 自己实现牛顿-柯特斯公式,求积公式的余项表达式,复合梯形公式,复合辛普森公式,龙贝格求积公式,自适应积分方法，高斯-勒让德求积公式,高斯-切比雪夫求积公式,多重积分
@Created By  : PyCharm
"""
from sympy.abc import x
from sympy import integrate, exp, diff
import numpy as np
from math import factorial


def trapz(a, b, intergrand):
    """
    实现梯形求积公式
    :param a: 区间左端点即起点a
    :param b: 区间右端点即终点b
    :param intergrand: 被积函数f(x)
    :return: 梯形求积积分
    """
    return (b - a) / 2 * (intergrand.subs(x, a) + intergrand.subs(x, b))


def newtown_cotes_integrate(equal_segment, interval_start, interval_end, symbol_t, f_x):
    """
    实现牛顿-柯特斯(Newton-Cotes)插值型积分
    :param symbol_t: 引进变换x=a+t*h后的积分变量
    :param equal_segment: 区间等分数n
    :param interval_start: 区间左端点即起点
    :param interval_end: 区间右端点即终点
    :param f_x: 被积函数 intergrand
    :return:牛顿-柯特斯(Newton-Cotes)插值积分
    """
    step_h = (interval_end - interval_start) / equal_segment
    f_k_array = np.array([f_x.subs(x, interval_start + k * step_h) for k in range(equal_segment + 1)])
    return (interval_end - interval_start) * np.sum(cotes_coefficient(equal_segment, symbol_t) * f_k_array)


def cotes_coefficient(equal_segment, symbol_t):
    """
    实现牛顿-柯特斯(Newton-Cotes)插值型求积公式的柯特斯系数
    :param equal_segment: 区间等分数n,其中1<=n<=7,n>7的牛顿-柯特斯公式计算不稳定，不使用.
    :param symbol_t: 引进变换x=a+t*h后的积分变量t
    :return: 柯特斯系数数组
    """
    if equal_segment not in range(1, 8):
        raise ValueError("Cotes coefficient must be an integer between 1 and 7")
    c_k_array = np.zeros(equal_segment + 1, dtype=np.float64)
    for k in range(equal_segment + 1):
        index = list(range(equal_segment + 1))
        index.pop(k)
        intergrand = np.prod(symbol_t - np.array(index))
        integral = integrate(intergrand, (symbol_t, 0, equal_segment))
        c_k_array[k] = (-1) ** (equal_segment - k) / (
                equal_segment * factorial(k) * factorial(equal_segment - k)) * integral
    return c_k_array


def quadrature_reminder(a, b, equal_segment, m, intergrand, a_k: np.ndarray):
    """
    实现[a,b]代数精度为m的求积公式余项表达式
    :param a: 区间左端点即起点a
    :param b: 区间右端点即终点b
    :param equal_segment: 区间等分数n
    :param m: 代数精度m
    :param intergrand: 被积函数f(x)
    :param a_k: 权系数Ak数组
    :return: 求积公式余项表达式
    """
    step_h = (b - a) / equal_segment
    f_x = x ** (m + 1)
    f_x_k_array = np.array([f_x.subs(x, a + k * step_h) for k in range(equal_segment + 1)])
    return 1 / factorial(m + 1) * (1 / (m + 2) * (b ** (m + 2) - a ** (m + 2)) - np.sum(a_k * f_x_k_array)) * diff(
        intergrand, x, m + 1), a, b


def com_trapz(a, b, equal_segment, intergrand):
    """
    实现复合梯形求积公式
    :param a: 区间左端点即起点a
    :param b: 区间右端点即终点b
    :param equal_segment: 区间等分数n
    :param intergrand: 被积函数f(x)
    :return: 复合梯形求积积分
    """
    step_h = (b - a) / equal_segment
    f_x_k_array = np.array([intergrand.subs(x, a + k * step_h) for k in range(1, equal_segment)])
    return step_h / 2 * (intergrand.subs(x, a) + 2 * np.sum(f_x_k_array) + intergrand.subs(x, b))


def com_simpson(a, b, equal_segment, intergrand):
    """
    实现复合辛普森求积公式
    :param a: 区间左端点即起点a
    :param b: 区间右端点即终点b
    :param equal_segment: 区间等分数n
    :param intergrand: 被积函数f(x)
    :return: 复合辛普森求积积分
    """
    step_h = (b - a) / equal_segment
    f_x_k1_array = np.array([intergrand.subs(x, a + (k+0.5) * step_h) for k in range(equal_segment)],
                            dtype=np.float64)
    f_x_k2_array = np.array([intergrand.subs(x, a + k * step_h) for k in range(1, equal_segment)],
                            dtype=np.float64)
    return step_h / 6 * (
            intergrand.subs(x, a) + 4 * np.sum(f_x_k1_array) + 2 * np.sum(f_x_k2_array) + intergrand.subs(x, b))


if __name__ == '__main__':
    # 辛普森公式及其余项表达式测试成功，来源详见来源详见李庆扬数值分析第5版P135,e.g.4
    """
    # 0.632333680003663
    inter_grand = exp(-x)
    print("辛普森公式积分为:{}".format(newtown_cotes_integrate(2, 0, 1, x, inter_grand)))
    # (-0.000347222222222221*exp(-x), 0, 1)
    print("辛普森公式为:{}".format(quadrature_reminder(a=0, b=1, equal_segment=2, m=3, intergrand=exp(-x), a_k=cotes_coefficient(2, x))[0]))
    # print(diff(exp(-x), x, 4) * (-1 / (180 * 16)))
    """
    # 复合梯形公式和复合辛普森公式测试成功，来源详见来源详见李庆扬数值分析第5版P135,e.g.2(1)
    """ 
    f_x = x / (4 + x ** 2)
    print("复合梯形公式积分:{}".format(com_trapz(a=0, b=1, equal_segment=8, intergrand=f_x)))  # 0.111402354529548
    print("复合辛普森公式积分:{}".format(com_simpson(a=0, b=1, equal_segment=8, intergrand=f_x)))  # 0.111571813252631
    """
